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ABSTRACT 

Regions of stellar and planetary interiors that are unstable according to the Schwarzschild criterion, 
but stable according to the Ledoux criterion, are subject to a form of oscillatory double- diffusive 
(ODD) convection often called "semi-convection" . In this series of papers, we use an extensive suite 
of three-dimensional (3D) numerical simulations to quantify the transport of heat and composition 
by ODD convection, and ultimately propose a new ID prescription that can be used in stellar and 
planetary structure and evolution models. The first paper in this series demonstrated that under cer- 
tain conditions ODD convection spontaneously transitions from an initially homogeneously turbulent 
state into a staircase of convective layers, which results in a substantial increase in the transport of 
heat and composition. Here, we present simulations of ODD convection in this layered regime, we 
describe the dynamical behavior of the layers, and we derive empirical scaling laws for the transport 
through layered convection. 
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I. INTRODUCTION 

A long-standing problem in modeling stellar and planetary interiors is how to quantify the transport of heat and com- 
position in regions that have a destabilizing (i.e., superadiabatic) thermal stratification and a s tabilizing compositional 
strat ification, such that the total density stratification is stabl e according to Ledoux's criterion (ISchwarzschild fc Harm 



1985 



1982; 



Leconte fc Chabrier 



19581 ). Such regions are present, for instance, in g iant planets (|Stevensonl l 

sive stars both main-sequ ence ( Merryfieldl 1995 ) and post- main-sequence ([Robertson fc Faulkner 



20121 ) , and in mas- 



Jliad; 



Langer et al 



Wooslev et al.ll2002l ). It has long been recognized that, because the molecular diffusivity of composition is typi- 



cally only a tiny fraction of the total (r adiative plu s conductive) diffusivity of temperatur e, these regions are subject 



to a form of double-diffusive instability ([Stern 



196C; 



Walin 



1964 : 



Veronis 



1965 



Kate 



1966f ). The turbulent convection 



that arises from this instability, which is often called "diffusive" or "oscillatory" convection, enhances the transport of 
heat and composition, changing the internal structure of the host object. 

Oscillatory double-diffusive convection is closely related to the phenomenon known as "semi-convection," but the 
latter is often also associated with the effects of convective overshoot and the dependence of thermal diffusivity 
on composition, neither of which are included in this study. Instead, we consider a localized region away from 

Itsw25® soe.ucsc.edul 



any radiative-convective boundary, in which the thermal and compositional diffusivities can be approximated as 
constant, and we study t he instability and turb ulence arising from the competition between thermal and compositional 
stratification. Following iPaparella et alj ( 20021 ) we refer to convection in this situation as "oscillatory double-diffusive 
(ODD) convection." 

The linear instability mechanism that leads to ODD convection can be understood from physical arguments. In 
the absence of thermal and compositional diffusion, vertically displaced parcels of fiuid oscillate adiabatically about 
their equilibrium position with a frequency determined by the total density stratification. But the presence of di ffusive 
proce sses can amplify the oscillation, producing an instability that resembles an overstable internal gravity wave (jKato 
19661 ). Indeed, if the thermal diffusivity is sufficiently large in comparison with the compositional diffusivity, then a 
displaced parcel will acquire heat (and therefore entropy) while below its equilibrium position, and thus will be slightly 
more buoyant than its surroundings when it returns to its equilibrium position. The amplitude of the oscillation grows 
exponentially until nonlinear effects are sufficient to saturate the instability, leading to a state of weakly turbulent 
ODD convection. 

After saturation, ODD convection initially takes the form of a homogenously turb ulent state, in whic h the time 
averaged stratification remains uniform. Previous work by 



Rosenblum et al 



( 2011 ) and 



Mirouh et al 



shown that the turbulen t transport of heat and composition in this homogeneous state is fairly weak. 



(2012) has 



Mirouh et al 



( 20121 hereafter 



Paper IT ) derived empirical formulae for the respective fluxes from an extensive suite of numerical 



simulations performed in local Cartesian domains. However, in certain cases the homogeneous state was found to 
spontaneously tra nsition i nto a new state of layered convection, for which the transport of heat and composition is 
greatly enhanced. iPaper 11 derived a simple semi-analytical criterion that predicts under which conditions spontaneous 
layering occurs. Under the physical conditions relevant to planetary and stellar interiors, spontaneous layering is 
expected within a significant fraction o f the dou ble- diffusive instability range, specifically in regions that are closer to 
being Ledoux- unstable (see Figure 1 in iPaper J ) . 

In this paper, we present further results from numerical simulations and use them to derive empirical formulae for 
the fluxes of heat and composition in layered convection. We determine how the fluxes depend on the thickness of 
the layers and on other governing parameters. We also describe some of the physical properties of layered convection, 
paying particular attention to the interfaces between the layers, and we compare our results with previous studies of 
ODD convection. Finally, we discuss our results in the light of recent layered models of giant planets. 



2. THE MATHEMATICAL MODEL AND NUMERICAL SCHEME 



The numerical model that we use to study ODD convection was described in detail in iPaper J . so here we mention 
only its most relevant features. The computational domain is a Cartesian box with uniform gravity g = —gez- We 
solve the evolution equations for momentum, temperature, and composition in the Boussinesq approximation, which 
are 
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'dt 



+ u • Vu 



— Vp ge^ 

Po Po 



dT 
'dt 



U-VT+ (To, - T^^)w = ktV^T, 



dt 



U • V/i + flQ^W = K^V^^, 



V-u = 0. 



(1) 
(2) 
(3) 
(4) 



Here p, p, T, and fi are perturbations of pressure, density, temperature, and mean molecular weight respectively, and 
w is the vertical component of the velocity u; all of these quantities are assumed to be periodic in all three Cartesian 
directions. We use subscript 'O's to denote the unperturbed, background fields, which we assume are functions of 
z alone. The background gradients of temperature, Tg^, and mean molecular weight, poz, are both negative in the 
situation considered here, as is the adiabatic temperature gradient. 



POz 



(5) 



ad 



We consider a small volume within which the diffusivities of momentum, temperature, kt, and composition, k^, 
can be treated as constant, and the density perturbations can be expressed as a linear function of the temperature 
and composition perturbations. 



P_ 

Po 



= l3fi — oT, 



where 



and /? 



1 

Po \dp 



(6) 
(7) 



Both a and (3 are positive for any realistic equation of state. We approximate a, (3, and the other coefficients 
in Equations (Il])-(l3]) as constant within our small domain, and we suppose that the background temperature and 
molecular weight profiles are linear in z, so that the total temperature and composition fields are Ttot — T + TqzZ and 
Mtot = + PozZ, with Tqz, Tq^, and poz all constant. The Schwarzschild and Ledoux stability criteria are then 



and 



a{Toz 
aiToz 



Tn^^) > 



Tof ) > ppoz 



(8) 
(9) 



respectively. We are concerned here with regions which satisfy ([9]) but not ([8]), i.e., regions in which Rq^ > 1, where 
R^^ is the "inverse density ratio," 



PIPQz 



(10) 



From a linear stability analysis of Equations (IT])-(l3]), it can be shown that such regions are unstable to double-diffusive 
instability if Rq^ lies in the range 



1 < i?Q ^ < -Rj. 



(Walin 



1964 : 



Kato 



19661) where R^ ^ is the critical value 



R^ — 



1 + Pr 
T + Pr' 



Here, Pr and r are the diffusivity ratios 



and 



Pr 



V j Kt 
i^^l/kt- 



(11) 
(12) 

(13) 
(14) 



In astrophysical objects, both Pr and r are typically of order 10 ^ or smaller, and so the instability window given by 
Equation (|11[) is rather wide. The instability is oscillatory, and has a typical horizontal lengthscale 

1/4 



ktv 



ag\Toz 



T^ad I 



(15) 



(e-g-, 



Baines h Gill 



19691) . 



We are interested primarily in the fluxes of heat and composition through ODD convection. As in iPaper I 
measure these fluxes in terms of Nusselt numbers, defined as 

(wT) 



we 



NU7 



1 



and 



Nu^ = 1 



(16) 
(17) 



where the notation ^ • ^ represents a volume average over the entire domain. The total vertical fluxes of heat, J-V, and 
composition, J-^, can be reconstructed from these as follows: 



PoCp 



= {wT) - ktTqz 
= -mrHTTo, + (Nut - l)«Trof 



and 



= -Nu^K^/io^, 



(18) 
(19) 



where Cp is the specific heat. We note that, in general, the thermal Nusselt number Nut defined by Equation (|16p 
is not the ratio of the total heat flux to the diff usive (i.e ., radiative plus conductive) heat flux, unless the adiabatic 
temperature gradient Tg^ is zero. However, as in Ipaper"^ . we can introduce the total "potential temperature" 



Aot =T+ (To, - )z 



(20) 



which has background gradient -iJoz = (Toz ^ ^^of ) and perturbation = T, in order to write the thermal Nusselt 
number as 



Nut = 1 



(21) 



So we may regard Nut as the ratio of the total and diffusive fluxes of potential temperature. Since and T are equal, 
in what follows we will refer to as simply "the temperature perturbation" . 
In double- diffusive convection, an important quantity is the "buoyancy flux ratio" 7~^, which we define as 



7-^ = 



gaKTt^OzNuT 



-iNu^ 



( Radko 



2003 



" Nut 

3. This is closely related to the core erosion "efficiency factor" x introduced by 



(22) 



Guillot et al 



( 2004 ) 



which measures the fraction of a planet's luminosity that is used to erode the stable compositional stratification at the 
boundary with the core. In the limit where the transport of both heat and composition is dominated by convection (i.e., 
Nut.^ 3> 1) this fraction is precisely 7^^. The gravitational potential energy released from the thermal stratification 
m ust excee d the potential energy used to erode the compositional stratification, which implies that 7^^ < 1. As shown 
in iPaper ll . 7"^ plays an importan t role in the formation of layers, and possibly also controls the evolution of layers 



after their formation ( Radko 



2Q05h 



For computational convenience, we work mainly with nondimensional quantities. As in IPaper 11 we nondimensionalize 

1 Unlike Radko 1 20031) . here we define 7 as the ratio of the total fluxes, rather than the turbulent fluxes. See lPaper 11 for more detail. 



Equations (HJ-® using the lengthscale d from Equation (IT5|) as well as the scales 

[t] = (f/KT, 



ad I 
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The governing equations are then 



T; 



■ad 
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where all variables are now dimensionless. The definitions p^ -([T7 | of the Nusselt numbers then become 

Nut = 1 



and 



Nu^ = 1 



We will also make use of alternative definitions of the Nusselt numbers, 

nuT = l + (|Vi9p) 

nu„ = 1 



and 



(23) 
(24) 
(25) 
(26) 

(27) 

(28) 
(29) 
(30) 

(31) 
(32) 

(33) 
(34) 



which are based on the rates of thermal and compositional dissipation, and are sometimes called "Cox numbers" . 
We note that uut and nu^ are guaranteed to be positive at all times, whereas Nut and Nu^ are not. However, it 
can be readily verified, from Equations (|28|) - (|30)) . that th e definitions (l3T]Hl32]) and (l33l)-(l34l) y ield the same result 

Shraiman &: Siggialll990l) . The advantage of 



1954 : 



after taking a time average in a statistically steady state (jMalkus 
retaining both forms is explained in Section [4.1l 

3. TYPICAL PROPERTIES OF LAYERED CONVECTION 

In principle, any particular i nstance of ODD convection is characterized by the three dimensionless parameters 
Rq^, Pr, and t. However, in Paper I it was shown that, within the instability window given by Equation (|lip . 
there is ano ther crit ical value of Rq^., say i?L^, below which ODD convection spontaneously transitions into layered 



convection. 



Paper I presented empirical scaling laws for i?^ ^ , and for the transport of heat and composition prior to 



the development of layers. It was shown that, in the parameter regime relevant to astrophysical objects, for which 
Pr ~ r ^ 1, the critical value ^ Pr^^^^. This result implies that there is a significant range of density ratios Rq^ 
for which layers form spontaneously, and so the transport through layer ed convection must be c onsidered in models 
of astrophysical objects!^ A preliminary study of layered convection by 



Rosenblum et al 



(j201ll) demonstrated that 

the turbulent transport is strongly dependent on the layer height, H say, in addition to Pr, r, and Rq^- Their results 
suggest that Nut and Nu^ both follow a power law in H, with an exponent between 1 and 4/3. In what follows we 

^ This is in contrast to the "fing ering" regime of do uble-diffusive convection, for which spontaneous layer format ion is not expected 



at astrophysical parameter values I Traxler et al.l l201ll ) , unless a different mechanism of layer formation operates I Brown et al 
submitted). 



20ia 



present a more comprehensive study of the dependence of the turbulent fluxes on each of the parameters i?, Pr, r, 
and Rq^- 

In order to provide a clear description of the transport properties of layered convection, and the subtleties involved 
in measuring the thermal and c ompositio nal fluxes in the layered phase, we first revisit one particular simulation 



that was previously presen ted in iPaper j . This simulation has i?g ^ = 1.5, Pr ~ 0.03, r — 0.03, and a domain size 
of (lOOd)'^. As described in lPaper j . this simulation spontaneously transitions from homogeneous ODD convection to 



layered convection, and the three layers that initially form subsequently merge into two layers, then into one layer. 
This evolution is illustrated in Figure [TJ which shows snapshots of the composition perturbation /i during each phase 
of the simulation, as well as profiles of the horizontally-averaged total "potential density," 



Ptot = Po(/3Mtot - a^^tot)- 



(35) 



The profiles show that ptot is roughly constant within each layer, but increases slightly with height. The local Ledoux 

d 

stability criterion can be expressed as ° < (cf. Equation (IS])), so this indicates that the layers are subject to 

oz 

overturning convection. Within each layer, the total composition /itot a-nd the total potential temperature iJtot are 
both roughly uniform, and there are strong jumps in /Ztot and ?9tot across the interfaces between adjacent layers. 




Ptot 




Ptot 



ptot 



ptot 



Figure 1. Composition perturbation fi (top row) and potential density profiles ptot(z) (bottom row) at four different times in a simulation 
with R. 







1.5, Pr = r = 0.03, and a domain size of (lOOd)^. The first column shows the homogeneous phase of ODD convection prior 
to layer formation. The remaining columns show layered convection with 3, 2, and 1 layers respectively. Because the perturbation is 
plotted in the top row, rather than the total composition, the layers appear as regions with a relatively weak vertical gradient, and the 
interfaces appear as regions with a strong vertical gradient of the opposite sign. The color bar for each plot is scaled to the maximum of 
the perturbation /i, which increases as the layers merge. In the profiles, ptot has been rescaled such that the total change in ptot across the 
domain is unity. 

Spontaneous layer formation similar to that observed here has previously been reported in simulations of double- 
diffusive convection in the "fingeri ng" regime, for which the thermal stratification is stabilizing a nd the compositiona l 



Stellmach et al 



20111) . but only for Prandtl numbers Pr > 1 (jTraxler et al 



2011) 



stratification is destabilizing (e.g.. 
However, in fingering convection the initial height of the layers is typically much larger than the lengthscale of the 
linear double-diffusive instability, whereas here we find that layers initially form on a scale comparable to that of the 
linear instability. In the simulation shown in Figure [l] for example, the fastest growing linear mode has a wavelength 



of roughly 20d, i.e., 1/5 of the domain size, and the initial layer height is 1/3 of the domain size. 

Figure m shows more detail of the interface structure during the 2-layer phase. We find that the interfaces are highly 
dynamic and turbulent, exhibiting structure on scales that are much smaller than the wavelength of the fastest growing 
linear mode. This is in marked contrast to the flat, laminar interfaces found in simulations and lab oratory experiments 
of ODD convection at larger Prandtl number (e.g. 



Noguchi fc Niino 



2010at 



Carpenter et al 



20121 both of which have 



Pr ~ 7) . The composition field exhibits structure on finer scales than the temperature field, as expected because of the 
smaller diffusivity of composition. On larger scales, the temperature and composition fields display similar features, 
such that the temperature field resembles a coarse-grained version of the composition field. By contrast, the density 
field has even finer structure, since it depends on the difference between the composition and temperature fields. 





Figure 2. Perturbations of temperature (left), composition (center), and density (right) from the simulation shown in Figure[T] at t = 1600. 
Regions where the perturbations are small have been made transparent in order to show the interfaces in more detail. 

The fluxes through this simulation, measured in terms of Nusselt numbers Nut.^ and nuT,;^, are plotted in Figure [3] 
as a function of time, along with the mean-squared velocity, (|up). As mentioned earlie r, the mean value of eac h flux 
increases significantly when layers first form, and increases still further when layers merge (jRosenblum et al Iboill The 
mean kinetic energy also increases, roughly in proportion with the fluxes, indicating that the simulation becomes more 
turbulent as the layers become thicker. We note that the two definitions for the Nusselt numbers both yield consistent 
results for the mean flux, although the definition based on the instantaneous fluxes exhibits larger fluctuations in the 
homogen eous, pre-layer phas e of the simulation, during which the dynamics are dominated by wave-like, oscillatory 



motions ( Mirouh et al 



2012 ) 



Between layer merger events the mean fluxes and kinetic energy oscillate quasi-periodically between high and low 
values, with relative departures from the mean flux of order unity. The oscillations occur as the flow in the simulation 
cycles between states that are highly turbulent and states that are more quiescent. In the quiescent state, the interfaces 
are roughly horizontal, whereas in the turbulent state the interfaces have undulations on a horizontal scale comparable 
to the thickness of the layers. For illustration. Figure 2] shows the composition perturbation, and the horizontally- 
averaged total composition /itot, at two times during the one-layer phase. At the first time, t — 1980, the interface is in 
a quiescent state, and at the second, t — 2020, the interface is in a more turbulent state. These two times correspond 
to minima and maxima in the Nusselt numbers, as shown in Figure [31 

The quasi-periodic nature of the oscillations in the fiuxes suggests that they are caused by large-scale gravity-wave 
motions of the interfaces. For an interfacial gravity wave of horizontal scale H , the oscillation frequency is 



' Aptot g 
Po H 



(36) 
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Figure 3. The Nusselt numbers, defined by Equations II31I I- II34II . and tiie mean kinetic energy in the simulation shown in Figures [Tl and [2l 
The initial linear instability grows exponentially until t fa 570, then saturates, leading to homogenous ODD convection. Layers first form 
at t 1200; the duration of the 3 layer, 2 layer, and 1 layer phases are indicated at the top of the plot. The vertical lines indicate the two 
times shown in Figure |4] Using nuy rather than Nuy to measure the mean fluxes has the advantage of filtering out fluctuations in the 
fluxes on short timescales, while retaining the important long timescale information. 



(jRavleighlll883l ). where Aptot represents the jump in potential density across the interface. For a staircase of layers of 
thickness H we have Aptot ~ Hpo\l3fioz — aiJozl, and so the characteristic gravity- wave frequency is 



w VgWik^T^^^l = N, (37) 

where N is the buoyancy frequency associated with the mean stratification, which in dimensionless units is iV = 
'\Jpt{Rq^ — 1). In the simulation shown in Figure |3] this corresponds to a (dimensionless) oscillation period of 27r/aj ~ 
50, which is approximately half of the typical period of the observed flux oscillations. We emphasize, however, that the 
disturbances to the interfaces seen in Figure HI are far more complex in structure than linear gravity waves. Moreover, 
the large-scale motions of the interfaces are not strictly periodic in time — they grow and decay cyclically, but they 
do not oscillate about zero. 



Previous models of layered convection in astrophysical objects fe.g.. ISpruit 



1992 ) have assumed that the interfaces 



are horizontal, laminar, and diffusive. If that were the case, the fluxes though the interfaces could be determined 
simply by measuring the vertical gradients of the horizontally-averaged Ttot and /itot profiles in the center of each 
interface. Figures [3] and HI show that this type of model cannot be valid here, because the fluxes are maximal when 
these gradients are weakest. In fact, since the interfaces are not horizontal, the fluxes depend on the local gradients 
within the interfaces and on the effective surface area of the interfaces. Figure S] suggests that variations in the surface 
area are the main cause of the oscillations in the fluxes. 
The interface behavior seen here somewhat resembles that seen in the compressible, two-dimensional simulations 
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t = 1980 




t = 2020 




-1.5 -1.0 -0.5 0.0 
/^tot 

Figure 4. The interface structure at times t = 1980 and t = 2020 during the one-layer phase. The profiles show the horizontal average of 
the total composition /itot, and the volume rendering shows the perturbation /i using the same color bar at both times. 



of iBielld (j200lf ) at low Prandtl number. iBielld also attributed the transport across the interfaces to the breaking 
of interfacial gravity waves, although the interpretation was complicated by interactions with the boundaries of the 
domain. The disturbances to the interface in Figure 0] have the same horizontal lengthscale as the width of the domain, 
lOOd, suggesting that, here also, the dynamics may be sensitive to the choice of domain size and the periodic boundary 
conditions. However, as described in Appendix |X1 similar behavior is seen in all of our layered simulations, and the 
results do not depend significantly on the size of the domain. 

4. MEAN FLUXES THROUGH LAYERED CONVECTION 

4.1. Method for Measuring the Mean Fluxes 

Our aim is to determine how the mean fluxes of heat and composition, measured in terms of dimensionless Nusselt 
numbers, depend on the dimensionless parameters Rq^, Pr, r, and the average (dimensionless) layer height, which we 
denote as H. As described in the last section, during each layered phase the fluxes of heat and composition oscillate 
quasi-periodically about a mean value, and reliably calculating the mean fluxes requires time averaging over several of 
these oscillations. We measure the mean fluxes and estimate the errors in our measurements as follows: 

1. We first identify the number of layers during each phase of the simulation, by counting the number of "steps" in 
the vertical profiles of the horizontally averaged fields, and matching each merger event with the corresponding 
increase in the Nusselt numbers. 

2. Once we have identified the time interval corresponding to a particular layered phase, we then divide this interval 
into four equal subintervals, and calculate the average values of Nut, Nu^, nu^, and nu^, during each subinterval 
using Equations piP " P4p . Thus we obtain eight separate measurements of the thermal and compositional fiuxes 
during each layered phase. 

3. We treat each of our eight measurements as independent, and calculate their mean and standard deviation, 
which yields estimates fo r the me an flux and the measurement error. (A similar method was used to estimate 



the measurement error in 



Paper j .) 



4. In the same fashion, we calculate the buoyancy flux ratio 7^^ for each quarter interval using Equation (|22p . and 
for both definitions of the Nusselt numbers. As before, we use the resulting eight measurements to estimate the 
mean value of 7^^ and our measurement error. 



10 

We have applied this procedure to each of the simulations from iPaper 11 in which layers formed spontaneously, as 
well as to many add itional simulations that have since been performed. The parameters of each simulation, including 
those from lPaper 11 that formed into layers, are listed in Table [T] 



We present below the results of these measurements and describe their dependence on the parameters Pr, r, Rq^ , 
and H . In principle these four parameters are sufficient to characterize any particular regime of layered convection. 
However, in our numerical simulations we also specify the size of the computational domain and the numerical resolu- 
tion, and so wc must check whether these additional parameters have any effect on the results. Appendix El describes 
some of the tests that we have performed to ensure that our measurements are robust. For reasons explained in the 
Appendix, only the flux measurements from simulations performed in d omains o f width ^ lOOd are presented below. 



We also note that the spontaneous emergence of layers described in iPaper 11 is not the only mechanism by which 



layers can develop in ODD convection. Laboratory experimen ts and observations of Ea r th's oceans indicat e that the 



formation of layers can be "triggered" by external influences (|Turner fc Stommell 11964 ; 



Kellev et al 



convection m ight even be metastable under conditions for which the linear ODD instability is absent (jVeronis 



20031). Lavered 



1965 



SpiegeJ[l972 ). In this paper, however, we only consider the regime in which layers form spontaneously, 1 < i?g < i?L^, 
in order to guarantee the reproducibility of our results. 

4.2. Parameter Dependence of the Mean Fluxes 

To facilitate comparison of our results with those of previous studies, we introduce the thermal Rayleigh number for 
layered convection, which is 



KTl' 



(38) 



where Hd is the (dimensional) layer height. Owing to computational limitations, our results cover only a limited range 
of Pr and t values. Moreover, by choosing to analyz e only th ose simulations in which layers form spontaneously, we 
are restricted to a rather narrow range for R^"^ (see IPaper I ). On the other hand, the range of Rar values that can 
be achieved in each simulation is limited only by the height of the computational domain. Moreover, the fact that 
layers merge allows us to test the dependence of the Nusselt numbers on Ra^ using data from a single simulation. The 
tallest simulation presented here has a domain of height 400d, and initially formed eight layers. Figure [5] shows the 
Nusselt numbers in this simulation, plotted as functions of time and Ra^. We find that both Nut and Nu^ are well 
approximated by a power law in Ra^ with an exponent of 1/3. This result can be explained physically if we suppose 
that the mean flux depends only on the jumps in potential temperature and composition across the interfaces, Ai9 and 
A/1, say, and not on any properties of the layers. More precisely, let us suppose that the (dimensional) buoyancy fluxes 
gaKT'&oz^'^T and (//SK^/iozNu^ depend only on the buoyancy jumps yaAi? and aPAn and th e diffusivities kt, k^, and 



On dimensional grounds, the fluxes must then be proportional to (Ai?)'*/^ ( Turner 



1965ft. In our simulations, the 



total change in temperature between the top and bottom of the domain is fixed, and so the jump across each interface 
is inversely proportional to the number of layers, and therefore proportional to the height of the layers. The exponent 
1/3 then follows from equation (1551) . 

We now ask whether all of our Nusselt number measurements can be fitted to the same power law. In order to fit 
all of the data simultaneously, we must make some assumption regarding the dependence not only on Ray, but also 
on Pr, T, and Rq^- We consider the thermal and compositional Nusselt numbers in turn, and then the flux ratio 7^^. 



4.2.1. The Heat Flux (Nut) 
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Table 1 

Parameters of each simulation 



Pr 


T 


i?o ^ 


Domain size 


Resolution 
(Fourier modes) 


Initial layer No. 


0.1 


0.1 


1.1 


100^ 


128^ 


2 


0.1 


0.1 


1.1 


100^ 


643 


2 


0.1 


0.1 


1.1 


50^ 


643 


1 


0.3 


0.3 


1.1 


10(f 


12S3 


2 


0.3 


0.3 


1.1 


100^ 


643 


2 


0.3 


0.3 


1.15 


10 ff^ 


64^ X 128 


2 


0.3 


0.3 


1.15 


100^ 


643 


2 


0.3 


0.3 


1.15 


200^x100 


1283 


2 


0.3 
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Figure 5. a) Time series of the Nusselt numbers from a simulation with Pr = t = 0.33, -Rg 



RaT 

1.15, and domain size 200^ X 400. The 



number of layers in each layered phase is indicated on the plot. The horizontal lines indicate the approximate duration of each layered 
phase and the mean flux determined by the method described in Section |4. II b) Dependence of Nut on Ray in the same simulation. The 
errors in the measurements were determined by the method described in Section |4. II The solid and dashed lines both have a slope of 1/3. 

Many analytical and heuristic arguments have been made that the turbulent transport of heat thr ough layered 



convection follows a p ower law in both the Rayleigh number, Ray, and the Prandtl number, Pr (e.g., ISpruit 



1992 



Balmforth et al. 



2006 ). But its dependence on the parameters ^ and r is less clear. As a first step toward fitting 

(39) 



our data to an analytical formula for Nut, we therefore assume a dependence of the form 

Nut-1 = Ra?,PrV(i?o\T), 



where the exponents a and 6, and the function /, are to be determined. We group all of our Nut measurements accord- 
ing to their values oiK^^ and r, which leads to 15 distinct groups: (i?o^^, r) = (1.1, 0.03), (1.1, 0.1), (1.1, 0.3), (1.15, 0.3), 
(1.15,0.33), (1.2,0.1), (1.2, 0.3), (1.25, 0.03), (1.25,0.1), (1.25,0.3), (1.4,0.1), (1.5,0.01), (1.5,0.03), (1.5,0.1), (2,0.01). 
We then look for a fit of the form 

NuT-1 = Ra?,PrVn, (40) 

where n = 1, 15 represents the index of each group. To determine the coefficients a, 6, and /„ we perform a weighted 
least-squares minimization with weights proportional to the inverse square of the estimated error in the measurements 
of Nut. The best fit to the data has exponents a — 0.34 ± 0.01 and h — 0.34 ± 0.03. The best fit values for the 
coefficients /„ are shown in Figure [51 along with a comparison between the data and the fit. Most of the data points 
lie within one standard error of the best fit line, although there are significant outliers for both low and high values 
of Nut. The lowest values of Nut are found in simulations that are only weakly turbulent, for which we would not 
expect a simple scaling law to hold. The presence outliers at the highest values of Nut, on the other hand, may result 
from insufficient time-averaging or numerical resolution in the simulations with the tallest layers. 

The fact that the best fit has exponents a h suggests that the heat flux depends only on the param e ters R n^ , 
r, and the product RaTPr. In that case the heat flux is independent of viscosity, as predicted bv ISpruitI (|l992l) for 
layered convection with RaT ^ 1 and Pr < 1. However, the model of ISpruitI (jl992r ) predicts exponents a = 6 = 1/4, 
rather than 1/3, which is incompatible with our results. 

1/3 

Using the result that Nut — 1 oc Pr ' we can directly compare simulations performed with different values of Pr 
but the same values of r and R^ ^ . An example is given in Figure [71 which overlays time series of Nut in three 
simulations that each have r = 0.3, R^^ — 1.2, and domain size (lOOd)'^. For these parameters, the best fit prediction 
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Figure 6. Left panel: The best-fit values for the coefficients /„ defined in Equation 1 140 I I. Right panel; Comparison between the measured 
values of Nuy — 1 and the best fit for the formula in Equation Ij40p v^ith exponents a = 0.34 and 6 = 0.34. In the online version of this 
figure, data from different parameter pairs (Rq^,t) are plotted in different colors. 

is Nut — 1 = 0.09i/^/"^Pr^/'^. Since each of these simulations has the same (dimensionless) domain height, we can 
directly compare their two-layer {H = 50) and one-layer {H — 100) phases; we find that, as predicted, (Nut — 1)/Pr^/'^ 
has the same dependence on H in each of these simulations. 

90 

80 

70 

" 60 
50 
40 
30 
20 
10 


3000 

Figure 7. Time series of the turbulent heat ffux, rescaled by Pr~^/^, in three simulations writh t = 0.3, Rq^ = 1.2, and domain height 
lOOd, but different values of Pr. The dashed lines indicate the flux predicted by the best fit for the two-layer and one-layer phases, 0.09H^^^ . 

The best- fit values for the coefficients /„, presented in Figure [SI show a weak but well-defined dependence on Rq^ 
and T (for fixed Ray and Pr). In particular, we find that decreasing either the composition gradient fioz or the 
compositional diffusivity produces a slight increase in the mean heat fiux. We interpret this result as evidence that 
changing fiQz and affects the heat fiux only indirectly, by changing the overall level of turbulence. 
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4.2.2. The Composttion Flux (IN^u^jj 
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We apply the same analysis for the compositional flux, i.e., we suppose that 

Nu^-1 = RajPrVn, 



(41) 



and fit the data to determine a, b, and /„. This time the best-fit values for the exponents are a — 0.37 ± 0.01 and 
b = 0.27 ± 0.04. This is roughly consistent with the trend observed in Figure[5]3, but indicates that the compositional 
flux, unlike the heat flux, depends on viscosity even for Pr ^ 1. This is perhaps not surprising given that the 
composition field is more sensitive than the temperature field to the small-scale motions of the fluid (e.g., see Figure [5]). 

Figure [5] shows the best fit values for the coefficients /„ and a comparison between the measured values of Nu^ 
and the best fit for Equation (|41l) . As found previously for Nut, there is a weak increase in Nu^ as Rq^ is made 
smaller, i.e., as the composition gradient is reduced. Unlike Nut^, however, Nu^ depends strongly on r, i.e., on the 
compositional diffusivity. In fact, we find that Nu^ is roughly inversely proportional to r, indicating that the mean 
turbulent flux of composition (wfi) is approximately independent of compositional diffusivity in the parameter regime 
considered here. 
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Figure 8. Left panel: The best-fit values for the coefficients /„ defined in Equation Hil l. The color bar has been rescaled to accommodate 
the larger range of values here. Right panel: Comparison between the measured values of Nu^ — 1 and the best fit formula in Equation l|41|l 
with exponents a = 0.37 and b = 0.27. 



4.2.3. The Buoyancy Flux Ratio (j ^) 

Taken at face value, the empirical power laws for Nut and Nu^ derived above imply that the flux ratio 7"^ oc 
Nu^/Nut oc Ra^"'^. However, as discussed in Section [21 7^^ cannot exceed unity in a statistically steady state, and 
therefore cannot follow a simple power law in Ray. As illustrated in Figure[9l we flnd that 7^^ is generally an increasing 
function of Ray, and possibly asymptotes to a constant value at large values of Rar- Further simulations at high Kelt 
are required to conflrm this trend, and to determine how this asymptotic value depends on the other parameters Rq^ , 
Pr, and r. Figure |9] suggests that 7"^ is a decreasing function of Rq^ , but again, more data is required to conflrm this 
suggestion. 

If 7~^ remains a monotonically increasing function of Ra^ for large Ra-r then there are important implications for 
the merging of layers, as we discuss in the next section. 



5. SUMMARY AND DISCUSSION 
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Figure 9. The dependence of on Ray for all simulations witti Pr = t = 0.3, and for tiie large-domain simulation shown in Figure [5l 
which has Pr = r = 0.33. Different colors and symbols indicate different values of Rq^ ■ 



We have measured the mean fluxes of heat and composition in an extensive suite of numerical simulations of layered 
ODD convection. In these simulations, layers develop spontaneously from an initially homogeneously turbulent state, 
and subsequently merge into a single layer spanning the entire height of the simulation domain. The formation and 
merging of layers is, to some extent, a stochastic process, but nevertheless occurs at a similar rate in simulations with 
different domain sizes and numerical resolutions. 

The development of layers leads to a significant increase in the mean fluxes, and the fluxes increase still further each 
time layers merge. We find that t he mean fluxe s are approximately proportional to the average thickness of the layers 
to the power 4/3, as predicted bv iTurneii (|l965l ). The heat flux is independent of viscosity, but the composition flux is 
not. Within the range of parameter values considered here, the thermal and compositional Nusselt numbers are well 
approximated by 

Nut - 1 = APr^/^Rai/^ 



and 



Nu„ 



1 



Sr-ipri/^Rai/' 



(42) 
(43) 



where the prefactors A and B are weakly dependent on the parameters Rq^ , r, Pr, and Ra^, and have typical values 
A ~ 0.1 and B ~ 0.03 respectively. The mean heat and composition flux can be calculated from these expressions 
using Equations (|18p and (jl9p . We emphasize, however, that the parameter values used in our simulations are rather 
far from true astrophysical values, and so any application of the formulae (|42p ~ (l43|) to astrophysical conditions requires 
an extrapolation beyond the range in which these formulae have been derived. 

In all our simulations the interfaces between the convective layers exhibit complex dynamical behavior that is very 
different from the laminar "interface scouring" and entra inment seen in studies of layered convection at higher Prandtl 
number (e.g., 



Linden &: Shirtcliffe 



1978 : 



Fernando 



19891 ). The structure of the interfaces, and their effective surface 



area, is highly time dependent, which produces order unity variations in the fluxes of heat and composition during 
each layered phase. Previous models for the transport by ODD convection have assumed that the interfaces remain 
flat and laminar, implying that the fluxes are limited by the vertical gradients of the horizontally averaged profiles 
within the interfaces. Our numerical results demonstrate that transport at low Prandtl number is more efficient than 
predicted by such quasi-one-dimensional models. We argue that the interface dynamics must be considered in future 
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models for transport under astrophysical parameter conditions. 

Since the transport by layered convection is strongly dependent on the layer thickness, we m ust consider what 
determines this thickness in stellar and planetary interiors. Recently, iLeconte fc Chabrieii (|2012| ) have proposed a 
model for the transport by ODD convection in Jupiter's interior. They use heuristic arguments to obtain upper 
and lower bounds for the height of layers in Jupiter, and thereby derive bounds for the thermal and compositional 
transport. In all of our simulations, however, the layers eventually merge into a single layer with the same height as the 
computational domain, so our results make no prediction for the height of layers in real objects. Whether simulations 
performed in ta ller domains w ould yield a saturation thickness for layers is unknown. For fingering convection, it has 
been argued by iRadkd (|2005r ) that there exists such a saturation thickness. The essence of his argument is that the 
buoyancy flux ratio 7"^ increases every time layers merge, but that there exists a maximum permitted value of 7"^, 
which is attained at a finite layer thickness H ~ iJo, say. Layers of thickness H ^ Hq therefore do not undergo further 
merging. In our simulations of ODD convection, by contrast, 7^^ appears to asymp tote toward a limiting val ue as 



Leconte &: Chabrier 



2012 ). An 



iJ — >■ 00, so this argument does not apply. (This is in contra st to a statement made by 
alternative picture for layer merging has been proposed bv iNoguchi fc Niind ( 2010bl) . but their model also predicts 
that the merging of layers proceeds indefinitely. 

We note that the Boussinesq approximation assumed in our numerical model breaks down once the layer thickness 
becomes comparable to the local pressure scale height of the planet or star. In practice, this may be what limits the 
size of layers in real astrophysical objects. At this point, convection within the layers would be almost indistinguishable 
from "normal" overturning convection, allowing the usual mixing-length approximation for the transport to be applied. 

An important issue that has not been explicitly addressed so far is the timescale for layer merging. Our results suggest 
that this ti mescale i s an i ncreasing function of layer height, which is in accordance with the theory of layer mergers 
proposed bv iRadkd (|2005l ). If the timescale for merging in an astrophysical object eventually becomes comparable to 
the global evolution timescale, then the layers will never reach a saturation thickness, and the evolution of the layers 
must then be modeled as an explicitly time-dependent process in ID structure and evolution models. Future papers 
will study the dynamics of layer mergers in more detail, in order to determine what sets the merging timescale, and 
how it depends on layer height. 

We thank B. R. Ruddick, E. A. Spiegel, and M. G. Wells for useful comments and suggestions. P. G. and T. W. 
were supported by funding from the NSF (NSF-0933759). Part of the computations were performed on the UGSC 
Pleiades supercomputer, purchased with an NSF-MRI grant. Others used computer resources at the National Energy 
Research Scientific Computing Genter (NERSG), which is supported by the Office of Science of the US Department 
of Energy under contract DE-AG03-76SF00098. 



APPENDIX 

INFLUENCE OF DOMAIN SIZE AND NUMERICAL RESOLUTION 

In this appendix we investigate the infiuence of numerical resolution and domain size on the initial layer height, the 
timescale for layer formation and merging, the time-averaged fluxes, and the amplitude of oscillations in the fluxes. 



Tests of Numerical Resolution 

Because of the disparity between the diffusivities of temperature and composition, double-diffusive convection always 
features a wide range of dynamical length scales, making this a tough problem to model numerically. Furthermore, 
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as mentioned in Section [31 the emergence of layers is associate d with an increase in the overall level of turbulence, 
and increased spectral energy at all scales in a given simulation ([Stellmach et al.ll2011l ). Properly resolving the layered 
phase therefore requires even greater numerical resolution than resolving the homogenous phase. In order to check 
that the thermal and compositional fluxes in our simulations are not sensitive to the numerical resolution, we have 
run several simulations that have the same physical parameters and domain size, but different numbers of grid points 
in the vertical and horizontal directions. 

As an example, Figure [10] compares time series of Nut.^ and nuy,^ from four simulations performed with different 
numerical resolutions. All four were performed in domains of size (50^)"^ with the same physical parameters as the 
simulation presented in Section [31 but with 48'^, 64'^, 96'^, and 128'^ Fourier modes respectively. The simulation 
with the coarsest resolution (48"^) only resolves structures on lengthscales > d, but even in this simulation the linear 
double-diffusive instability is well resolved (recall that the fastest growing linear mode has wavelength ~ 20c? at these 
parameters). Each simulation undergoes a homogeneous phase, then forms three layers, which later merge into two 
layers, then one layer. The duration of each layered phase varies between the four simulations, suggesting that the 
merging of layers is, to some extent, a stochastic process. However, the four simulations produce similar fluxes during 
each phase, and the variations in each flux about its mean are comparable in all four cases. 
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Figure 10. Comparison of Nussclt numbers between simulations with 48 , 
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96 , and 128 Fourier modes in each spatial direction. 



All four cases have Pr : 



: 0.03, Rn 



1.5, and a domain size of (SOd)^. The two rows show the Nusselt numbers calculated from the 



instantaneous fluxes II31I I and II32I I. and the instantaneous dissipation rates II33I I and 1134 l l. 



In each simulation, the most stringent test of the numerical resolution occurs during the final merging event, which 
marks the start of the one-layer phase. Figure [TT] compares the density fields at this point in the lowest and highest 
resolution simulations shown in Figure [TOl The lowest resolution simulation is clearly under-resolved, with significant 
features at the Nyquist scale, whereas the highest resolution simulation displays no numerical artifacts. Nevertheless, 
the mean fluxes through these two simulations are approximately equal in the one-layer phase, which suggests that 
the transport of both heat and composition is dominated by scales > d in these simulations. 
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Figure 11. Horizontal cross-sections of density p from the 48'' resolution simulation at t = 1530 (left) and the 128'' resolution simulation 
at t = 1800 (right). Both cases show cross-sections through the interface at the start of the one-layer phase, which is typically the least 
well-resolved stage in any simulation. 

Tests of Domain Size and Aspect Ratio 

Figure [T^] compares time series of the Nusselt numbers from four simulations performed in domains of different 
heights and widths. Each simulation has Pr = r = 0.3, -Rq"^ = 1.15, and the same numerical resolution, i.e., the same 
(dimensionless) distance between neighboring grid points. These four simulations exhibit the following features, which 
we find to be typical: 

• The cases with wider domains (200c?) take longer to form layers than the cases with narrower domains (lOOd). 

• The initial layer height is not very sensitive to the size of the domain, provided that the domain is sufficiently 
wide. In Figure \\% for e xample, the initia l layer height is SOd in all four cases. However, in other simulations 



(not shown here, but see 



Rosenblum et al 



(j2011l) for instance) it has been found that in very narrow domains 
layers form much earlier, and the initial layer height is smaller. This result indicates that the horizontally- 
periodic boundary conditions, which encourage horizontal coherence in the simulation, can artificially accelerate 
the formation of layers when the domain is sufficiently narrow. For this reason, we only use data from simulations 
that were performed in domains of width ^ lOOd when calculating the mean fluxes in Section W7\\ 

• The mean fluxes of heat and composition, for layers of a particular height, are consistent between runs with 
different domain sizes. 

• During each layered phase, the cases with wider domains exhibit smaller variations in the value of each flux 
about its mean value. This suggests that the variations in the fluxes are horizontally localized, and therefore 
have less effect on the mean flux in wider domains. Similarly, the simulation in the tallest domain (of height 
400d) exhibits the smallest variations in the fluxes during each layered phase, because in this simulation the 
volumetric mean flux represents an average over a greater number of layers and interfaces. 

Figure [T^ also illustrates a practical difficulty that arises when comparing results from simulations performed in 
domains of different heights, which is that simulations in taller domains usually undergo more layered phases than 
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Figure 12. Comparison of Nusselt numbers between simulations in domains of different sizes. All four cases have Pr = r = 0.3 and 



simulations in shorter domains, and the length of each layered phase is correspondingly shorter. For example, the 
simulation with a domain height of 400c? passes through phases in which the average layer height is 400(i/7 and 400d/6; 
these states cannot be achieved in the simulations with a domain height of lOOd. As a result, the time series of the 
Nusselt numbers in the taller simulation are less obviously step-like, and the mean flux during each layered phase is 
more difficult to measure accurately. This also makes it difficult to determine whether the domain height affects the 
rate at which layers merge. The full time series of both Nusselt numbers from this simulation are shown in Figure [5^. 
We find that the time between layer merger e vents increase s as the number of layers decreases, which is in accordance 
with the model of layer merging proposed by iRadkd (j2005[ ). 

Another difficulty associated with measuring fluxes in simulations with multiple layers is that the layers are not 
necessarily all of the same height. In that case, a complete description of the state of the system would require 
(at least) measuring the height of each layer, and the jumps in temperature and composition across each interface, 
potentially introducing many new parameters on which the mean flux might depend. Furthermore, in such cases we 
would not expect the system to be in statistical equilibrium, and so the mean flux might not be well deflned. To 
avoid dealing with this additional complexity, in this study we assume that, in each layered phase, all the layers have 
approximately the same height, while recognizing that this assumption introduces an additional source of error into 
our flux measurements. 
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